%g gives optimal lending policy, index keeps track of the index of optimal
%lending in the mspace vector, p is optimal prices
%K - 1 if propping-up effect is active. 0 if not. 
%N - number of large lenders
%alpha - vector of share of non-borrowers in good and bad state
%pl_share - vector of private-label share in good and bad state
%alpha_h - vector of high-quality borrower share in good and bad state
%phi_bh - vector of high-quality borrower endowment shocks in good and bad state
%phi_bl - vector of low-quality borrower endowment shocks in good and bad state
%q_s - transition matrix [HH HL; LH LL]
%beta - beta of representative bank (beta_hat in paper)
%delta - housing depreciation 
%c - housing supply parameter
%e - borrower endowments when positive


function [g,index,p] =  vfi_benchmark(K,N,alpha,pl_share, alpha_h, phi_bh, phi_bl, q_s, beta,delta,c,e)
%K is 1 if propping-up 
%K is 0 is not active
alpha_l = 1-pl_share-alpha-alpha_h; %vector of share of low-income borrowers in high and low state
outside_d(1) = alpha(1)+pl_share(1); %vector of outside housing demand in high state
outside_d(2) = alpha(2)+pl_share(2); %vector of outside housing demand in low state

%Grid for m_ts 
mspace=linspace(0,1,100); %rowvector of M_t
%%
n=100;
% For each state, I am building a matrix which gives the value for every
% feasible m and m^{- combination}
u=cell(2,1);
for l=1:2 %loops over two states
    u{l}=-Inf*ones(n,n);
%create a transition matrix for the utility of each feasible m' for every m
for i=1:n %loop over m
    for j=1:n %loop over m'
        P(j) = c*(mspace(j)+outside_d(l));%Price if bank makes m' loans
  if K==1 %propping-up active 
    u{l}(i,j)= mspace(i)*(1-delta)*P(j)-(mspace(j)*P(j))-c*(min(mspace(j),alpha_h(l)))^2/2*(N-1)-c*(min(alpha_l(l),max(mspace(j)-alpha_h(l),0)))^2/2*(N-1)-c*mspace(j)*(N-1)*outside_d(l)+beta*min(mspace(j),alpha_h(l))*e*(q_s(l,1)*phi_bh(1) + q_s(l,2)*phi_bh(2))+beta*max(mspace(j)-alpha_h(l),0)*e*(q_s(l,1)*phi_bl(1) + q_s(l,2)*phi_bl(2)) -(N-1)*min(mspace(j),alpha_h(l))*(min(alpha_l(l),max(mspace(j)-alpha_h(l),0)))*c;   % end
  else %propping-up inactive
    u{l}(i,j)= -(mspace(j)*P(j))-c*(min(mspace(j),alpha_h(l)))^2/2*(N-1)-c*(min(alpha_l(l),max(mspace(j)-alpha_h(l),0)))^2/2*(N-1)-c*mspace(j)*(N-1)*outside_d(l)+beta*min(mspace(j),alpha_h(l))*e*(q_s(l,1)*phi_bh(1) + q_s(l,2)*phi_bh(2))+beta*max(mspace(j)-alpha_h(l),0)*e*(q_s(l,1)*phi_bl(1) + q_s(l,2)*phi_bl(2)) -(N-1)*min(mspace(j),alpha_h(l))*(min(alpha_l(l),max(mspace(j)-alpha_h(l),0)))*c;   % end
  end
     if mspace(j)>alpha_h(l)+alpha_l(l)
        u{l}(i,j) = -Inf;
     end
    end
end
end

v=zeros(n,2); %create value function with initial guess 0 for all Is
tv=ones(n,2);
index=ones(n,2); %index keeps track of the index for optimal lending
G=ones(n,2); %G gives optimal lending
tolerance=1*10^(-1);
count=0;
while max(abs(tv-v))>tolerance
      if count==0
          v=zeros(n,2);
      else
          v=tv;
      end
for i=1:2 %loops over states
    for k=1:n %loop over m outstanding loans
        p_temp(k,1) = c*(outside_d(1)+G(k,1));
        p_temp(k,2) = c*(outside_d(2)+G(k,2));
        %propping up
        if K==1 %when propping-up is active, u (above) includes price so that bank accounts for its effect on price
        E_p(k,i) = 0; 
       %no propping-up
        else %when propping-up is inactive, u (above) does not include price. This implies that the bank is taking expected future prices as given and not internalizing its effect on them.
        E_p(k,i) = N* min(mspace(k),alpha_h(i)+alpha_l(i)) *(q_s(i,1)*p_temp(k,1)+q_s(i,2)*p_temp(k,2)); 
        end
        [tv(k,i) ,index(k,i)]=max(u{i}(k,:)+beta/N*(1-delta)*E_p(k,i)+beta*q_s(i,:)*v'); 
         G(k,i)=mspace(index(k,i));
         g(k,i)= mspace(index(k,i))/N;
    end
    count=count+1;
end

    
end
p(:,1) = c*(outside_d(1)+G(:,1)); %p is optimal prices
p(:,2) = c*(outside_d(2)+G(:,2));
end